/*
 * NativeLapack.cs
 *
 * Wrapper around dnA.Native.ILapack.
 *
 * Copyright (c) 2005, dnAnalytics. All rights reserved.
 */

#region Using Directives
using System;
using dnAnalytics;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Native;
#endregion Using Directives

namespace dnAnalytics.LinearAlgebra
{
    internal sealed class NativeLapack : ISimplifiedLapack
    {
        #region Fields
        private ILapack lapack = Provider.LapackInstance;
        private IBlas blas = Provider.BlasInstance;
        #endregion Fields

        #region Public Methods
        public void LUFactor(int order, double[] factor, int[] pivots)
        {
            int[] ipiv;
            lapack.Dgetrf(order, order, factor, order, out ipiv);
            for (int i = 0; i < pivots.Length; i++)
            {
                pivots[i] = ipiv[i] - 1;
            }
        }

        public void LUInverse(int order, double[] inverse, int[] pivots)
        {
            int[] ipiv = new int[order];
            Array.Copy(pivots, ipiv, pivots.Length);
            for (int i = 0; i < ipiv.Length; ++i)
            {
                ipiv[i] += 1;
            }

            lapack.Dgetri(order, inverse, order, ipiv);
        }

        public void LUSolve(int order, int columns, double[] factor, int[] pivots, double[] rightSide)
        {
            int[] ipiv = new int[order];
            Array.Copy(pivots, ipiv, pivots.Length);
            for (int i = 0; i < ipiv.Length; ++i)
            {
                ipiv[i] += 1;
            }
            lapack.Dgetrs('N', order, columns, factor, order, ipiv, rightSide, order);
        }

        public double Norm(char norm, int rows, int columns, double[] data)
        {
            return lapack.Dlange(norm, rows, columns, data, rows);
        }

        public double ConditionNumber(char norm, int order, double anorm, double[] luData, int[] pivots)
        {
            double rcond = 0;
            lapack.Dgecon(norm, order, luData, order, anorm, ref rcond);
            return rcond;
        }

        public int CholeskyFactor(int order, double[] data)
        {
            int ret = lapack.Dpotrf('L', order, data, order);
            for (int i = 0; i < order; i++)
            {
                for (int j = 0; j < order && i > j; j++)
                {
                    data[i * order + j] = 0;
                }
            }
            return ret;
        }

        public void CholeskySolve(int order, int columns, double[] data, double[] rhs)
        {
            lapack.Dpotrs('L', order, columns, data, order, rhs, order);
        }

        //uses dgeqrf, we should switch to dgeqp3.
        public void QRFactor(int m, int n, double[] q, double[] r, double[] tau)
        {
            lapack.Dgeqrf(m, n, r, m, tau);

            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < m && j < n; j++)
                {
                    if (i > j)
                    {
                        q[j * m + i] = r[j * m + i];
                    }
                }
            }

            //compute the q elements explicitly
            if (m <= n)
            {
                lapack.Dorgqr(m, m, m, q, m, tau);
            }
            else
            {
                lapack.Dorgqr(m, m, n, q, m, tau);
            }
        }

        public void QRSolve(int m, int n, int bn, double[] q, double[] r, double[] c, double[] tau)
        {
            lapack.Dormqr('L', 'T', m, bn, n, r, m, tau, c, m);
            blas.Dtrsm(BlasOrderType.Column, BlasSideType.Left, BlasUpLoType.Upper, BlasTransType.NoTrans,
                 BlasDiagType.NonUnit,
                 n, bn, 1, r, m, c, m);
        }

        public int SVDDecomposition(int m, int n, double[] a, double[] s, double[] u, double[] v)
        {
            char jobu = 'A';
            char jobvt = 'A';
            if (u == null)
            {
                jobu = 'N';
                jobvt = 'N';
            }
            return lapack.Dgesvd(jobu, jobvt, m, n, a, m, s, u, m, v, n);
        }

        #endregion Public Methods
    }
}